##################################################
# R Script for: A Replication of 'Representative Bureaucracy and the Willingness to Coproduce'

# Please note: This R Script is taken from the long-term archived CodeOcean Capsule of the research project. 
# You can easily reproduce the results if you access the script from the capsule: 

# https://codeocean.com/capsule/2129531/tree/v1

# The CodeOcean Capsule contains the used R Environment and likely saves you some hassle with the R packages and library commands. 


# 1. Setup ####################################################################

# 1.1 Packages ================================================================

library(dplyr)
library(stargazer)
library(Hmisc)
library(ggpubr)
library(grid)
library(gridExtra)
library(TOSTER)
library(AER)
library(table1)
library(htmltools)
library(MASS)

# 1.2 Functions ===============================================================
#IVRegress

ivregress = function(second, first, data){
  s.names = all.vars(second)
  f.names  = all.vars(first)
  data.names = names(data)
  all.names = c(s.names,f.names)
  resp   = s.names[1]
  endog  = f.names[1]
  inst   = f.names[-1]
  explan = s.names[-1]
  exog   = explan[explan!=endog]
  exog.f = paste(exog,collapse="+")
  inst.f = paste(inst, collapse="+")
  RHS    = paste(exog.f, inst.f, sep="+")
  first.form = as.formula( paste(endog, "~", RHS))
  first.lm = lm(first.form, data)
  ftest    = linearHypothesis(first.lm, inst, rep(0,length(inst)))
  x.hat    = fitted(first.lm)
  data2    = cbind(data,x.hat)
  iname    = paste(endog,".hat",sep="")
  names(data2) = c(names(data), iname)
  data2.names = names(data2)
  RHS2 = paste(exog.f,iname,sep="+")
  second.form = as.formula(paste(resp, "~", RHS2))
  second.lm = lm(second.form, data2)
  Xmat  = data2[c(exog,endog)]
  Xmat2  = data2[c(exog,iname)]
  z = summary(second.lm)
  X     = as.matrix(cbind(1,Xmat))
  X2    = as.matrix(cbind(1,Xmat2))
  Y     = data[resp]
  fit   = X%*%second.lm$coefficients
  res   = Y - fit
  xPx     = t(X2)%*%X2
  xPx.inv = solve(xPx)
  z$cov.unscaled     = xPx.inv
  z$residuals        = res
  z$sigma            = sqrt(mean(res^2))
  varcovmat          = z$cov.unscaled*z$sigma
  coef               = z$coefficients[,1]
  IV.SE              = z$sigma*sqrt(diag(xPx.inv))
  t.iv               = coef/IV.SE
  p.val              = 2*(1-pnorm(abs(t.iv)))
  z$coefficients[,2] = IV.SE
  z$coefficients[,3] = t.iv
  z$coefficients[,4] = p.val
  result = list(summary(first.lm),z,ftest)
  names(result) = c("first", "second", "ftest")
  print("IV object successfully created. Use sum.iv() on object")
  print("to learn about your 2SLS Regression")
  return(invisible(result))
  }
sum.iv = function(reg.iv, first=FALSE,
                  second=TRUE, ftest=FALSE) {
    x= rep(0,2)
    if(first==TRUE) x[1] = 1
    if(second==TRUE) x[2]= 2
    if(ftest==TRUE) x[3]= 3
    print(reg.iv[x])
  }

# Function that estimates the conditinal effects for certain values of the moderator
conditional.effects <-  function(model, x, m, quantiles = 3){
  interact = paste0(x,':',m)
  beta.hat = coef(model) 
  covs = vcov(model)
  z0  = quantile(model$model[,m], seq(0 , 1, 1/quantiles))
  dy.dx = beta.hat[x] + beta.hat[interact]*z0
  se.dy.dx = sqrt(covs[x, x] + z0^2*covs[interact, interact] + 2*z0*covs[x, interact])
  upr = dy.dx+1.96*se.dy.dx
  lwr = dy.dx-1.96*se.dy.dx
  data.frame(m = z0, b = dy.dx, lwr, upr)
}

# Build mean indices from multiple variables
mean_index <- function (df, name, vars) {
  M1 <- dplyr::select(df, vars)
  M2 <- rowMeans(M1, na.rm = TRUE)
  M2 <- tibble::tibble(M2)
  colnames(M2) <- name
  df <- dplyr::bind_cols(df, M2)
  return(df)
}

# Calculate standard error
se_func <- function(var) {
  sd <- sd(var)
  n <- length(var)
  se <- sd / sqrt(n)
  return(se)
}

# Calculate confidence interval (lower bound)
lower_ci_func <- function(var) {
  m <- mean(var)
  sd <- sd(var)
  n <- length(var)
  se <- sd / sqrt(n)
  lower_ci <- m - qt(1 - (0.05 / 2), n - 1) * se
  return(lower_ci)
}

# Calculate confidence interval (upper bound)
upper_ci_func <- function(var) {
  m <- mean(var)
  sd <- sd(var)
  n <- length(var)
  se <- sd / sqrt(n)
  lower_ci <- m + qt(1 - (0.05 / 2), n - 1) * se
  return(lower_ci)
}

# format p values in APA format
pformat <- function(p) {
  paste0("p ",
         ifelse(p >= 0.001,
                paste0("= ",
                       weights::rd(p, 3)),
                paste0("< .001")))
}

# Chi2 results in APA format
chi_result <- function (x, y) 
{
  suppressWarnings(chi_result <- stats::chisq.test(x, y))
  p <- pformat(chi_result[["p.value"]]) 
  result <- 
    paste0("Chi2(", 
           as.integer(chi_result[["parameter"]][["df"]]),
           ") = ",
           broman::myround(
             chi_result[["statistic"]][["X-squared"]], 
             1),
           ", ",
           p)
  return(result)
}

# ANOVA result in APA format
aov_result <- function (aov_object) {
  anova_result <- summary(aov_object)
  p <- pformat(anova_result[[1]][["Pr(>F)"]][1])
  result <- paste0("F(", as.integer(anova_result[[1]][["Df"]][1]), 
                   ",", as.integer(anova_result[[1]][["Df"]][2]), ") = ", 
                   broman::myround(anova_result[[1]][["F value"]][1], 2), 
                   ", ", p)
  return(result)
}

# Print all columns of a tibble
tibble_print_vars <- function (tibble) 
{
  default <- getOption("tibble.width")
  options(tibble.width = Inf)
  print(tibble)
  options(tibble.width = default)
}

data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}

# Correlation table
# x is a matrix containing the data
# method : correlation method. "pearson"" or "spearman"" is supported
# removeTriangle : remove upper or lower triangle# results :  if "html" or "latex"
# the results will be displayed in html or latex format# labels_rows and labels_cols are character vectors for labeling rows and columns
corstars <- function(x, method=c("pearson", "spearman"),
                     removeTriangle=c("upper", "lower"),
                     result=c("none", "html", "latex"),
                     labels_rows, labels_cols, 
                     cap = c("Correlation"), filename = ""){
  #Compute correlation matrix
  x <- as.matrix(x)
  correlation_matrix<-Hmisc::rcorr(x, type=method[1])
  R <- correlation_matrix$r # Matrix of correlation coeficients
  p <- correlation_matrix$P # Matrix of p-value 
  
  ## Define notions for significance levels; spacing is important.
  mystars <- ifelse(p < .01,"**", ifelse(p < .05, "* ", "  "))
  #mystars <- ifelse(p < .001, "****", ifelse(p < .001, "*** ", ifelse(p < .01, "**  ", ifelse(p < .05, "*   ", "    "))))
  
  ## trunctuate the correlation matrix to two decimal
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1]
  
  ## build a new matrix that includes the correlations with their apropriate stars
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
  diag(Rnew) <- paste(diag(R), " ", sep="")
  rownames(Rnew) <- colnames(x)
  colnames(Rnew) <- paste(colnames(x), "", sep="")
  
  ## remove upper triangle of correlation matrix
  if(removeTriangle[1]=="upper"){
    Rnew <- as.matrix(Rnew)
    Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
    Rnew <- as.data.frame(Rnew)
  }
  
  ## remove lower triangle of correlation matrix
  else if(removeTriangle[1]=="lower"){
    Rnew <- as.matrix(Rnew)
    Rnew[lower.tri(Rnew, diag = TRUE)] <- ""
    Rnew <- as.data.frame(Rnew)
  }
  
  ## remove last column and return the correlation matrix
  Rnew <- cbind(Rnew[1:length(Rnew)-1])
  rownames(Rnew) <- labels_rows
  colnames(Rnew) <- labels_cols
  if (result[1]=="none") return(Rnew)
  else{
    if(result[1]=="html") print(xtable::xtable(Rnew, caption = cap), 
                                type="html",
                                caption.placement = "top",
                                file = filename)
    else print(xtable::xtable(Rnew, caption = cap), 
               type="latex")
  }
} 

# 2. Import data ##############################################################

# 2.1 Dataset =================================================================

ds <- read.csv('../data/SymbolicRepresentation.csv')

# datafile denotes as -> ds

summary(ds)
ds_r <- ds

# 3. Data management ##########################################################

# 3.1 rename variables ========================================================

ds_represent <- ds_r %>% rename(rand_org = ds.RD02,
                                rand_exp = ds.RD03,
                                leg1 = ds.PM01_01,
                                leg2 = ds.PM01_02,
                                leg3 = ds.PM01_03,
                                leg4 = ds.PM01_04,
                                leg5 = ds.PM01_05,
                                leg6 = ds.PM01_06,
                                leg7 = ds.PM01_07,
                                leg8 = ds.PM01_08,
                                leg9 = ds.PM01_09,
                                leg10 = ds.PM01_10,
                                WLG1 = ds.DT01_01,
                                WLG2 = ds.DT01_02,
                                WLG3 = ds.DT01_03,
                                man = ds.DT05,
                                gender = ds.DE01,
                                politic = ds.DE04,
                                PSM1 = ds.DE07_01,
                                PSM2 = ds.DE07_02,
                                PSM3 = ds.DE07_03,
                                PSM4 = ds.DE07_04,
                                work = ds.DE10,
                                education = ds.DE11,
                                age = ds.DE14_01,
                                BAW = ds.DE15)


# 3.2 Add Dummies for Treatments ============================================

ds_represent$rand_exp[ds_represent$rand_exp == 1] <- 1
ds_represent$rand_exp[ds_represent$rand_exp == 2] <- 0
ds_represent$rand_exp[ds_represent$rand_exp == 3] <- 2

ds_represent$rand_exp <- as.factor(ds_represent$rand_exp)

# 3.3 Manipulation Check ====================================================


ds_represent$man_num <- as.numeric(as.factor(ds_represent$man))

ds_represent$man_num[ds_represent$man_num == 1] <- 0
ds_represent$man_num[ds_represent$man_num == 2] <- 1
ds_represent$man_num[ds_represent$man_num == 3] <- 2
ds_represent$man_num[ds_represent$man_num == 4] <- 3
ds_represent$man_num[ds_represent$man_num == 5] <- 4

mcheck <- lm(man_num ~ rand_exp, data = ds_represent)
print(mcheck)
summary(mcheck)

group_by(ds_represent, rand_exp) %>%
  summarise(
    count = n(),
    mean = mean(man_num, na.rm = TRUE),
    sd = sd(man_num, na.rm = TRUE)
  )

ggline(ds_represent, x = "rand_exp", y = "man_num", 
       add = c("mean_se", "jitter"), 
       order = c("0", "1", "2"),
       ylab = "Number of female officials", xlab = "Treatment group")

res.aov <- aov(man_num ~ as.factor(rand_exp), data = ds_represent)
summary(res.aov)

TukeyHSD(res.aov)

stargazer(mcheck, type = "html", digits = 2, out = "../results/manip.htm",
          star.cutoffs = c(0.05, 0.01, 0.001))

# 3.4 Legitimacy Scale ====================================================

# 3.4.1 Explorative Factor Analyses =======================================

psych::fa.parallel(ds_represent %>% dplyr::select(leg1,
                                           leg2,
                                           leg3,
                                           leg4,
                                           leg5,
                                           leg6,
                                           leg7,
                                           leg8,
                                           leg9,
                                           leg10),
                   fa = "fa")

fa.legit <- psych::fa(ds_represent %>%
                        dplyr::select(leg1,
                               leg2,
                               leg3,
                               leg4,
                               leg5,
                               leg6,
                               leg7,
                               leg8,
                               leg9,
                               leg10),
                      nfactors = 3, 
                      fm = "pa", rotate = "promax",
                      scores = "regression")
print(fa.legit, digits = 2, sort = TRUE)

print(fa.legit$loadings, cutoff=.3)

# 3.4.2 Build Indices =======================================

ds_represent <- mean_index(ds_represent, "legit_pragmatic", 
                      c("leg1",
                        "leg2",
                        "leg3"))

ds_represent <- mean_index(ds_represent, "legit_moral", 
                      c("leg4",
                        "leg5",
                        "leg6",
                        "leg7"))

ds_represent <- mean_index(ds_represent, "legit_cognitive", 
                      c("leg8",
                        "leg9",
                        "leg10"))

psych::alpha(ds_represent %>% dplyr::select(leg1,
                                     leg2, 
                                     leg3))

psych::alpha(ds_represent %>% dplyr::select(leg4,
                                     leg5,
                                     leg6,
                                     leg7))

psych::alpha(ds_represent %>% dplyr::select(leg8,
                                     leg9,
                                     leg10))

# 3.5 Political Orientation =======================================

ds_represent$politic <- as.numeric(as.factor(ds_represent$politic))

# 3.5 Public Service Motivation =======================================

# 3.5.1 Explorative Factor Analyses =======================================

psych::fa.parallel(ds_represent %>%
                     dplyr::select(PSM1,
                            PSM2,
                            PSM3,
                            PSM4))

fa.PSM <- psych::fa(ds_represent %>%
                      dplyr::select(PSM1,
                             PSM2,
                             PSM3,
                             PSM4),
                    nfactors = 1, 
                    fm = "minres", rotate = "oblimin",
                    scores = "regression")
print(fa.PSM, digits = 2, sort = TRUE)

psych::alpha(ds_represent %>% dplyr::select(PSM1,
                                PSM2,
                                PSM3,
                                PSM4))

# 3.5.2 Build Index ==============================================

ds_represent <- mean_index(ds_represent, "PSM_mean", 
                      c("PSM1",
                        "PSM2",
                        "PSM3",
                        "PSM4"))

# 4 Sample Description ##############################################################

# 4.1 Sample Characteristics (Using original wide data format) ======================

# Number of Participants
nrow(ds_represent)

# Gender (factor)
ds_represent <- ds_represent %>% mutate(gender = factor(gender, labels = c("male", "female")))

# Gender (Dummy 0/1)
ds_represent <- ds_represent %>% mutate(gender_dummy = as.numeric(gender) - 1)

# PSM
mean(ds_represent$PSM_mean, na.rm =TRUE)

# Age
mean(ds_represent$age, na.rm = TRUE)

ds_represent$age <- as.numeric(ds_represent$age)

ds_represent$agegroup <- cut(ds_represent$age, breaks=c(17, 29, 40, 50, 60, 70), right = FALSE)
ds_represent$agegroup

ds_represent$agegroup1 <- dplyr::between(ds_represent$age, 17, 28)
ds_represent$agegroup2 <- dplyr::between(ds_represent$age, 29, 39)
ds_represent$agegroup3 <- dplyr::between(ds_represent$age, 40, 49)
ds_represent$agegroup4 <- dplyr::between(ds_represent$age, 50, 59)
ds_represent$agegroup5 <- dplyr::between(ds_represent$age, 60, 70)

summary(ds_represent$agegroup)

# Work

ds_represent$worknumeric <- as.numeric(as.factor(ds_represent$work))

ds_represent$work_priv[ds_represent$worknumeric == 1] <- 1
ds_represent$work_priv[ds_represent$worknumeric != 1] <- 0

ds_represent$work_npo[ds_represent$worknumeric == 2] <- 1
ds_represent$work_npo[ds_represent$worknumeric != 2] <- 0

ds_represent$work_oef[ds_represent$worknumeric == 3] <- 1
ds_represent$work_oef[ds_represent$worknumeric != 3] <- 0

ds_represent$work_abl[ds_represent$worknumeric == 4] <- 1
ds_represent$work_abl[ds_represent$worknumeric != 4] <- 0

# Education

ds_represent$edunumeric <- as.numeric(as.factor(ds_represent$education))

ds_represent$edu_7[ds_represent$edunumeric == 1] <- 1
ds_represent$edu_7[ds_represent$edunumeric != 1] <- 0

ds_represent$edu_HSA[ds_represent$edunumeric == 2] <- 1
ds_represent$edu_HSA[ds_represent$edunumeric != 2] <- 0

ds_represent$edu_RSA[ds_represent$edunumeric == 3] <- 1
ds_represent$edu_RSA[ds_represent$edunumeric != 3] <- 0

ds_represent$edu_FHR[ds_represent$edunumeric == 4] <- 1
ds_represent$edu_FHR[ds_represent$edunumeric != 4] <- 0

ds_represent$edu_AB[ds_represent$edunumeric == 5] <- 1
ds_represent$edu_AB[ds_represent$edunumeric != 5] <- 0

# Gender
prop.table(table(ds$gender))

# 4.2 Sample Characteristics Table =========================================

# 4.3 Summary Table (Alternative Code) ========================================

table1::label(ds_represent$gender) <- "Gender"
table1::label(ds_represent$age) <- "Age"
table1::label(ds_represent$agegroup1) <- "18-29"
table1::label(ds_represent$agegroup2) <- "30-39"
table1::label(ds_represent$agegroup3) <- "40-49"
table1::label(ds_represent$agegroup4) <- "50-59"
table1::label(ds_represent$agegroup5) <- "60-69"
table1::label(ds_represent$work_priv) <- "Employment Private"
table1::label(ds_represent$work_oef) <- "Employment Public"
table1::label(ds_represent$work_npo) <- "Employment NPO"
table1::label(ds_represent$work_abl) <- "Unemployed"
table1::label(ds_represent$edu_7) <- "Less than 7 years"
table1::label(ds_represent$edu_HSA) <- "CSE"
table1::label(ds_represent$edu_RSA) <- "GCSE"
table1::label(ds_represent$edu_FHR) <- "High School (FH Reife)"
table1::label(ds_represent$edu_AB) <- "High School (Abitur)"
table1::label(ds_represent$politic) <- "Political Orientation"

table <- table1::table1(~ gender + age + agegroup1 + agegroup2 + agegroup3 +
                 agegroup4 + agegroup5 + work_priv + work_oef +
                 work_npo + work_abl + edu_7 + edu_HSA + edu_RSA +
                 edu_FHR + edu_AB + politic | rand_exp, data = ds_represent, 
               output = "html", export = "../results/sampleoverview")

save_html(table, file = "../results/sampleoverview.htm", background = "white", libdir = "lib", lang = "en")

# 5 Main Analysis ##############################################################

# 5.1 Regression Analysis ======================================================

overall1 <- lm(WLG1 ~ rand_exp, data = ds_represent)
summary(overall1)

overall2 <- lm(WLG2 ~ rand_exp, data = ds_represent)
summary(overall2)

overall3 <- lm(WLG3 ~ rand_exp, data = ds_represent)
summary(overall3)

overall_g <- lm(WLG1 ~ rand_exp + gender + politic + 
                  work + education + PSM_mean + 
                  age, data = ds_represent)

summary(overall_g)

overall_mod <- lm(WLG1 ~ rand_exp + gender + politic + 
                    work + education + PSM_mean + 
                    age + legit_pragmatic + legit_moral + 
                    legit_cognitive, data = ds_represent)

summary(overall_mod)

overall_mod2 <- lm(WLG1 ~ rand_exp + gender + 
                     politic + work + education + 
                     PSM_mean + age + legit_pragmatic + 
                     legit_moral + legit_cognitive +
                     rand_exp*legit_pragmatic + 
                     rand_exp*legit_moral + 
                     rand_exp*legit_cognitive,
                   data = ds_represent)

summary(overall_mod2)

# 5.2 Dependent Variables Cut-Off ======================================================

summary(ds_represent$WLG1)
summary(ds_represent$WLG2)
summary(ds_represent$WLG3)

ds_represent$WLG1_d[ds_represent$WLG1 > 3] <- 1
ds_represent$WLG1_d[ds_represent$WLG1 <= 3] <- 0

ds_represent$WLG2_d[ds_represent$WLG2 > 3] <- 1
ds_represent$WLG2_d[ds_represent$WLG2 <= 3] <- 0

ds_represent$WLG3_d[ds_represent$WLG3 > 3] <- 1
ds_represent$WLG3_d[ds_represent$WLG3 <= 3] <- 0

# 5.3 Sub-Sample Female ======================================================

ds_represent$gender <- as.numeric(ds_represent$gender)

ds_women <- ds_represent[ which(ds_represent$gender == 2), ]

# 5.4 Sub-Sample Male ======================================================

ds_male <- ds_represent[ which(ds_represent$gender == 1), ]

# 5.5 Repeated Regression Sub-Samples ======================================

overall_w <- lm(WLG1 ~ rand_exp, data = ds_women)
summary(overall_w)

overall_w2 <- lm(WLG2 ~ rand_exp, data = ds_women)
summary(overall_w2)

overall_w3 <- lm(WLG3 ~ rand_exp, data = ds_women)
summary(overall_w3)

overall_m <- lm(WLG1 ~ rand_exp, data = ds_male)
summary(overall_m)

overall_m2 <- lm(WLG2 ~ rand_exp, data = ds_male)
summary(overall_m2)

overall_m3 <- lm(WLG3 ~ rand_exp, data = ds_male)
summary(overall_m3)

# 5.6 Extended Regression (+ Controls) =========================================

overall_g_w <- lm(WLG1 ~ rand_exp + gender + 
                    politic + work + education + 
                    PSM_mean + age, data = ds_represent)

summary(overall_g_w)

overall_mod_w <- lm(WLG1 ~ rand_exp + gender + 
                      politic + work + education + 
                      PSM_mean + age + legit_pragmatic + 
                      legit_moral + legit_cognitive, data = ds_represent)
summary(overall_mod_w)

overall_mod2_w <- lm(WLG1 ~ rand_exp + gender + 
                       politic + work + education + 
                       PSM_mean + age + legit_pragmatic + 
                       legit_moral + legit_cognitive +
                       rand_exp*legit_pragmatic + 
                       rand_exp*legit_moral + 
                       rand_exp*legit_cognitive, 
                  data = ds_represent)

summary(overall_mod2_w)

overall_g_w2 <- lm(WLG2 ~ rand_exp + gender + 
                     politic + work + education + 
                     PSM_mean + age, data = ds_represent)

summary(overall_g_w2)

overall_mod_w2 <- lm(WLG2 ~ rand_exp + gender + 
                       politic + work + education + 
                       PSM_mean + age + legit_pragmatic + 
                       legit_moral + legit_cognitive, data = ds_represent)

summary(overall_mod_w2)

overall_mod2_w2 <- lm(WLG2 ~ rand_exp + gender + 
                        politic + work + education +
                        PSM_mean + age + legit_pragmatic + 
                        legit_moral + legit_cognitive +
                        rand_exp*legit_pragmatic + 
                        rand_exp*legit_moral + 
                        rand_exp*legit_cognitive, 
                     data = ds_represent)

summary(overall_mod2_w2)

overall_g_w3 <- lm(WLG3 ~ rand_exp + gender + 
                     politic + work + education + 
                     PSM_mean + age, data = ds_represent)

summary(overall_g_w3)

overall_mod_w3 <- lm(WLG3 ~ rand_exp + gender + 
                       politic + work + education + 
                       PSM_mean + age + legit_pragmatic + 
                       legit_moral + legit_cognitive, data = ds_represent)

summary(overall_mod_w3)

overall_mod2_w3 <- lm(WLG3 ~ rand_exp + gender + 
                        politic + work + education + 
                        PSM_mean + age + legit_pragmatic + 
                        legit_moral + legit_cognitive +
                        rand_exp*legit_pragmatic + 
                        rand_exp*legit_moral + 
                        rand_exp*legit_cognitive, 
                     data = ds_represent)

summary(overall_mod2_w3)

# 5.7 OLS Regression - Dichotomized =========================================

overall_d <- glm(WLG1_d ~ rand_exp, data = ds_represent)
summary(overall_d)

overall_d2 <- glm(WLG2_d ~ rand_exp, data = ds_represent)
summary(overall_d2)

overall_d3 <- glm(WLG3_d ~ rand_exp, data = ds_represent)
summary(overall_d3)

overall_g_d <- lm(WLG1_d ~ rand_exp + gender + 
                    politic + work + education + 
                    PSM_mean + age, data = ds_represent)

summary(overall_g_d)

overall_mod_d <- lm(WLG1_d ~ rand_exp + gender + 
                      politic + work + education +
                      PSM_mean + age + legit_pragmatic + 
                      legit_moral + legit_cognitive, data = ds_represent)

summary(overall_mod_d)

overall_mod2_d <- lm(WLG1_d ~ rand_exp + gender + 
                       politic + work + education +
                       PSM_mean + age + legit_pragmatic + 
                       legit_moral + legit_cognitive +
                       rand_exp*legit_pragmatic + 
                       rand_exp*legit_moral + 
                       rand_exp*legit_cognitive, 
                   data = ds_represent)

summary(overall_mod2_d)

# 5.8 OLS Regression - Dichotomized (Sub-Samples) =========================================

overall_w_d <- lm(WLG1_d ~ rand_exp, data = ds_women)
summary(overall_w_d)

overall_g_w_d <- lm(WLG1_d ~ rand_exp + gender + 
                      politic + work + education + 
                      PSM_mean + age, data = ds_represent)

summary(overall_g_w_d)

overall_mod_w <- lm(WLG1 ~ rand_exp + gender + 
                      politic + work + education + 
                      PSM_mean + age + legit_pragmatic + 
                      legit_moral + legit_cognitive, data = ds_represent)

summary(overall_mod_w)

overall_mod2_w <- lm(WLG1 ~ rand_exp + gender + 
                       politic + work + education + 
                       PSM_mean + age + legit_pragmatic + 
                       legit_moral + legit_cognitive +
                       rand_exp*legit_pragmatic + 
                       rand_exp*legit_moral + 
                       rand_exp*legit_cognitive, 
                     data = ds_represent)

summary(overall_mod2_w)

# 6 Graphics for Results ##############################################################

# 6.1 Randomization Checks ===========================================================

# Test-Statistics
aov_age <- aov(age ~ rand_exp, data = ds_represent)
aov_politic <- aov(politic ~ rand_exp, data = ds_represent)

ds_represent$work <- as.numeric(as.factor(ds_represent$work))
ds_represent$education <- as.numeric(as.factor(ds_represent$education))

aov_work <- aov(work ~ rand_exp, data = ds_represent)
aov_education <- aov(education ~ rand_exp, data = ds_represent)

chi_gender <- chiSquare(gender ~ rand_exp, data = ds_represent)

summary(aov_age)
summary(aov_politic)

summary(aov_work)
summary(aov_education)

TukeyHSD(aov_age)
TukeyHSD(aov_politic)
TukeyHSD(aov_work)
TukeyHSD(aov_education)

chi_gender

ds_represent$rand_exp_num <- as.numeric(ds_represent$rand_exp) 

# 6.2 Correlation Table ===========================================================

# Prepare data
df_cor <- ds_represent %>%
  dplyr::select(WLG1, WLG2, WLG3, age , gender, politic, work, education,
           PSM_mean, legit_pragmatic, legit_moral, legit_cognitive)
df_cor <- as.data.frame(df_cor)
df_cor$work <- as.numeric(df_cor$work)
df_cor$education <- as.numeric(df_cor$education)
df_cor$gender <- as.numeric(df_cor$gender)

# Table
corstars(df_cor, method="pearson", removeTriangle="upper", result="html",
         cap = "Correlations matrix", filename = "../results/corr.html",
         labels_rows = c("(1) Willingness (Lay Judge)", 
                         "(2) Willingness (Probationer)", 
                         "(3) Willingness (Correspondent)", 
                         "(4) Age",
                         "(5) Gender", 
                         "(6) Political Orientation", 
                         "(7) Employment Sector",
                         "(8) Educational Level",
                         "(9) Public Service Motivation",
                         "(10) Pragmatic Legitimacy",
                         "(11) Moral Legitimacy",
                         "(12) Cognitive Legitimacy"),
         labels_cols = 1:11)

# 6.3 Regression Analysis ===========================================================

# Regression Total Sample

stargazer(overall1, overall_d, overall2, overall_d2, overall3, overall_d3, 
          type = "html", digits = 2, summary = F, 
          dep.var.labels = c("Lay judge",
                             "Lay judge (>3)",
                             "Probationer",
                             "Probationer (>3)",
                             "Correspondent",
                             "Correspondent (>3)"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          omit.stat = "F", star.cutoffs = c(0.05, 0.01, 0.001),
          out ="../results/regression_table1.htm")

# Regression Gender Subgroups

stargazer(overall_w, overall_m, overall_w2, overall_m2, overall_w3, overall_m3, 
          type = "html", digits = 2, summary = F, 
          dep.var.labels = c("Lay judge",
                             "Probationer",
                             "Correspondent"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          omit.stat = "F", star.cutoffs = c(0.05, 0.01, 0.001),
          out ="../results/regression_table2.htm")

# 6.4 Visualization ===========================================================

## Lay Judge

ds_represent$rand_exp <- factor(x = ds_represent$rand_exp, labels = c("All male", "Female/male", "All female"))
ds_represent$gender <- factor(x = ds_represent$gender, labels = c("Men", "Women"))

ds_summary <- data_summary(ds_represent, varname = "WLG1", 
                        groupnames = c("rand_exp", "gender"))

j_ds <- ggplot(ds_summary, aes(x = rand_exp, y = WLG1, group = gender, color = gender)) +
  ylim(0, 5) +
  geom_line(aes(linetype = gender), position=position_dodge(0.3)) +
  geom_point(aes(shape = gender, color = gender), fill= "white", size = 6, position=position_dodge(0.3)) +
  geom_errorbar(aes(ymin = WLG1 - 1*sd, ymax = WLG1 + 1*sd),
                width=.2,
                position=position_dodge(0.3)) + 
  labs(x= "", y = "Willingness (Lay judge)") +
  theme_classic() + 
  theme(legend.title=element_blank(), legend.position = c(0.2, 0.9)) +
  scale_color_grey()

ggsave("../results/WLG1.png", j_ds, width = 4, dpi = 1000)

## Probationer

ds_summary2 <- data_summary(ds_represent, varname = "WLG2", 
                           groupnames = c("rand_exp", "gender"))

j_ds2 <- ggplot(ds_summary2, aes(x = rand_exp, y = WLG2, group = gender, color = gender)) +
  ylim(0, 5) +
  geom_line(aes(linetype = gender), position=position_dodge(0.3)) +
  geom_point(aes(shape = gender, color = gender), fill= "white", size = 6, position=position_dodge(0.3)) +
  geom_errorbar(aes(ymin = WLG2 - 1*sd, ymax = WLG2 + 1*sd),
                width=.2,
                position=position_dodge(0.3)) + 
  labs(x= "", y = "Willingness (Probationer)") +
  theme_classic() + 
  theme(legend.title=element_blank(), legend.position = "none") +
  scale_color_grey()

ggsave("../results/WLG2.png", j_ds2, width = 4, dpi = 1000)

## Correspondent

ds_summary3 <- data_summary(ds_represent, varname = "WLG3", 
                           groupnames = c("rand_exp", "gender"))

j_ds3 <- ggplot(ds_summary3, aes(x = rand_exp, y = WLG3, group = gender, color = gender)) +
  ylim(0, 5) +
  geom_point(aes(shape = gender, color = gender), fill= "white", size = 6, position=position_dodge(0.3)) +
  geom_line(aes(linetype = gender), position=position_dodge(0.3)) +
  geom_errorbar(aes(ymin = WLG3-1*sd, ymax = WLG3 + 1*sd),
                width=.2,
                position=position_dodge(0.3)) + 
  labs(x= "", y = "Willingness (Correspondent)") +
  theme_classic() + 
  theme(legend.title=element_blank(), legend.position = "none") +
  scale_color_grey()

ggsave("../results/WLG3.png", j_ds3, width = 4, dpi = 1000)

## Save as one Graphic

grid.arrange(j_ds, j_ds2, j_ds3, ncol = 3)

g <- arrangeGrob(j_ds, j_ds2, j_ds3, ncol = 3) #generates g
ggsave(file="../results/results.png", g, width = 9, dpi = 600)

# 6.5 Additional Analysis: TOST ===========================================================

psych::describeBy(ds_represent$WLG1, group = ds_represent$rand_exp)

WLG1a <- TOSTtwo(m1=2.66, m2=2.67, sd1=1.28, sd2=1.36, 
        n1=334, n2=329, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05)

WLG1b <- TOSTtwo(m1=2.66, m2=2.7, sd1=1.28, sd2=1.32, 
        n1=334, n2=337, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05)

psych::describeBy(ds_represent$WLG2, group = ds_represent$rand_exp)

WLG2a <- TOSTtwo(m1=2.49, m2=2.33, sd1=1.19, sd2=1.23, 
        n1=334, n2=329, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05)

WLG2b <- TOSTtwo(m1=2.49, m2=2.44, sd1=1.19, sd2=1.17, 
        n1=334, n2=337, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05)

psych::describeBy(ds_represent$WLG3, group = ds_represent$rand_exp)

WLG3a <- TOSTtwo(m1=2.42, m2=2.36, sd1=1.26, sd2=1.33, 
        n1=334, n2=329, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05)

WLG3b <- TOSTtwo(m1=2.42, m2=2.4, sd1=1.26, sd2=1.21, 
        n1=334, n2=337, low_eqbound_d = -0.2, 
        high_eqbound_d = 0.2, alpha = 0.05, plot = T)

table(WLG1a)

# 6.6 Robustness Checks ===========================================================

# Analysis as Treated - Sub Groups

ds_treat <- ds_represent

ds_treat$rand_exp <- as.numeric(ds_treat$rand_exp)

ds_treat$success[ds_treat$man_num == 0 & ds_treat$rand_exp == 1] <- 1
ds_treat$success[ds_treat$man_num == 2 & ds_treat$rand_exp == 2] <- 1 
ds_treat$success[ds_treat$man_num == 4 & ds_treat$rand_exp == 3] <- 1 

ds_treat$success[is.na(ds_treat$success)] <- 0  

summary(ds_treat$success)  
table(ds_treat$success)  

ds_treat$success_a <- 1

ds_treat$success_a[ds_treat$man_num >= 2 & ds_treat$rand_exp == 1] <- 0
ds_treat$success_a[ds_treat$man_num == 0 & ds_treat$rand_exp == 2] <- 0
ds_treat$success_a[ds_treat$man_num == 4 & ds_treat$rand_exp == 2] <- 0
ds_treat$success_a[ds_treat$man_num >= 3 & ds_treat$rand_exp == 3] <- 0

table(ds_treat$success_a)  

ds_treat_s <- ds_treat[ which(ds_treat$success == 1), ]
ds_treat_sa <- ds_treat[ which(ds_treat$success_a == 1), ]

ds_treat_s$rand <- as.factor(ds_treat_s$rand_exp)
ds_treat_sa$rand <- as.factor(ds_treat_sa$rand_exp)

### Regression Analysis 

# Main Effects (Conservative)
  
aat1 <- lm(WLG1 ~ rand, data = ds_treat_s)
summary(aat1)

aat2 <- lm(WLG2 ~ rand, data = ds_treat_s)
summary(aat2)

aat3 <- lm(WLG3 ~ rand, data = ds_treat_s)
summary(aat3)

aat_g <- lm(WLG1 ~ rand + gender + politic + 
                  work + education + PSM_mean + 
                  age, data = ds_treat_s)

summary(aat_g)

aat_mod <- lm(WLG1 ~ rand + gender + politic + 
                    work + education + PSM_mean + 
                    age + legit_pragmatic + legit_moral + 
                    legit_cognitive, data = ds_treat_s)

summary(aat_mod)

aat_mod2 <- lm(WLG1 ~ rand + gender + 
                     politic + work + education + 
                     PSM_mean + age + legit_pragmatic + 
                     legit_moral + legit_cognitive +
                     rand_exp*legit_pragmatic + 
                     rand_exp*legit_moral + 
                     rand_exp*legit_cognitive,
                   data = ds_treat_s)

summary(aat_mod2)

# Main Effects (Less Conservative)

aat1 <- lm(WLG1 ~ rand, data = ds_treat_sa)
summary(aat1)

aat2 <- lm(WLG2 ~ rand, data = ds_treat_sa)
summary(aat2)

aat3 <- lm(WLG3 ~ rand, data = ds_treat_sa)
summary(aat3)

aat_g <- lm(WLG1 ~ rand + gender + politic + 
              work + education + PSM_mean + 
              age, data = ds_treat_sa)

summary(aat_g)

aat_mod <- lm(WLG1 ~ rand + gender + politic + 
                work + education + PSM_mean + 
                age + legit_pragmatic + legit_moral + 
                legit_cognitive, data = ds_treat_sa)

summary(aat_mod)

aat_mod2 <- lm(WLG1 ~ rand + gender + 
                 politic + work + education + 
                 PSM_mean + age + legit_pragmatic + 
                 legit_moral + legit_cognitive +
                 rand_exp*legit_pragmatic + 
                 rand_exp*legit_moral + 
                 rand_exp*legit_cognitive,
               data = ds_treat_sa)

summary(aat_mod2)

# Gender Subgroups

ds_female_s <- ds_treat_sa[ which(ds_treat_sa$gender_dummy == 1),]
ds_male_s <- ds_treat_sa[ which(ds_treat_sa$gender_dummy == 0),]

# Main Effects Female & Male

f_aat1 <- lm(WLG1 ~ rand, data = ds_female_s)
summary(f_aat1)

f_aat2 <- lm(WLG2 ~ rand, data = ds_female_s)
summary(f_aat2)

f_aat3 <- lm(WLG3 ~ rand, data = ds_female_s)
summary(f_aat3)

m_aat1 <- lm(WLG1 ~ rand, data = ds_male_s)
summary(m_aat1)

m_aat2 <- lm(WLG2 ~ rand, data = ds_male_s)
summary(m_aat2)

m_aat3 <- lm(WLG3 ~ rand, data = ds_male_s)
summary(m_aat3)

# 6.7 CACE Analysis - 2SLS ===========================================================

# Separating Treatment

table(ds_treat$success_a)

ds_treat$mixed[ds_treat$rand_exp == 2] <- 1
ds_treat$mixed[ds_treat$rand_exp != 2] <- 0

ds_treat$female[ds_treat$rand_exp == 3] <- 1
ds_treat$female[ds_treat$rand_exp != 3] <- 0

ivmodel1 <- ivreg(formula = WLG1 ~ rand_exp | success_a, data = ds_treat)
ivmodel2 <- ivreg(formula = WLG2 ~ rand_exp | success_a, data = ds_treat)
ivmodel3 <- ivreg(formula = WLG3 ~ rand_exp | success_a, data = ds_treat)

summary(ivmodel1, vcov = sandwich, diagnostics = T)
summary(ivmodel1)

summary(ivmodel2, vcov = sandwich, diagnostics = T)
summary(ivmodel2)

summary(ivmodel3, vcov = sandwich, diagnostics = T)
summary(ivmodel3)

# 2SLS for Gender Sub Sample Female

ivmodel4 <- ivreg(formula = WLG1 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 1)
ivmodel5 <- ivreg(formula = WLG2 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 1)
ivmodel6 <- ivreg(formula = WLG3 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 1)

summary(ivmodel4, vcov = sandwich, diagnostics = T)
summary(ivmodel4)

summary(ivmodel5, vcov = sandwich, diagnostics = T)
summary(ivmodel5)

summary(ivmodel6, vcov = sandwich, diagnostics = T)
summary(ivmodel6)

# 2SLS for Gender Sub Sample Male

ivmodel7 <- ivreg(formula = WLG1 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 0)
ivmodel8 <- ivreg(formula = WLG2 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 0)
ivmodel9 <- ivreg(formula = WLG3 ~ rand_exp | success_a, data = ds_treat, subset = gender_dummy == 0)

summary(ivmodel7, vcov = sandwich, diagnostics = T)
summary(ivmodel7)

summary(ivmodel8, vcov = sandwich, diagnostics = T)
summary(ivmodel8)

summary(ivmodel9, vcov = sandwich, diagnostics = T)
summary(ivmodel9)

# 7. Reviewer Requests ###################################################

# 7.1 Ordered Logistic Regression ========================================

# 7.1.1 Full Sample ======================================================

ds_represent$WLG1_f <- as.factor(ds_represent$WLG1)
ds_represent$WLG2_f <- as.factor(ds_represent$WLG2)
ds_represent$WLG3_f <- as.factor(ds_represent$WLG3)

overall_olr_1 = polr(WLG1_f ~ rand_exp, data = ds_represent, Hess = TRUE)
summary(overall_olr_1)

overall_olr_2 = polr(WLG2_f ~ rand_exp, data = ds_represent, Hess = TRUE)
summary(overall_olr_2)

overall_olr_3 = polr(WLG3_f ~ rand_exp, data = ds_represent, Hess = TRUE)
summary(overall_olr_3)

summary_table1 <- coef(summary(overall_olr_1))
pval <- pnorm(abs(summary_table1[, "t value"]),lower.tail = FALSE)* 2
summary_table1 <- cbind(summary_table1, "p value" = round(pval,3))
summary_table1

summary_table2 <- coef(summary(overall_olr_2))
pval <- pnorm(abs(summary_table2[, "t value"]),lower.tail = FALSE)* 2
summary_table2 <- cbind(summary_table2, "p value" = round(pval,3))
summary_table2

summary_table3 <- coef(summary(overall_olr_3))
pval <- pnorm(abs(summary_table3[, "t value"]),lower.tail = FALSE)* 2
summary_table3 <- cbind(summary_table3, "p value" = round(pval,3))
summary_table3

overall_olr_1$AIC <- AIC(overall_olr_1)
overall_olr_2$AIC <- AIC(overall_olr_2)
overall_olr_3$AIC <- AIC(overall_olr_3)

stargazer(overall_olr_1, overall_olr_2, overall_olr_3, 
          type = "html", digits = 2, summary = F, 
          dep.var.labels = c("Lay judge",
                             "Probationer",
                             "Correspondent"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          omit.stat = "F", star.cutoffs = c(0.05, 0.01, 0.001),
          out ="../results/regression_robustness1.htm")

# 7.1.1 Sub Groups =======================================================

# Male

ds_male$WLG1_f <- as.factor(ds_male$WLG1)
ds_male$WLG2_f <- as.factor(ds_male$WLG2)
ds_male$WLG3_f <- as.factor(ds_male$WLG3)

overall_olr_1_m = polr(WLG1_f ~ rand_exp, data = ds_male, Hess = TRUE)
summary(overall_olr_1_m)

overall_olr_2_m = polr(WLG2_f ~ rand_exp, data = ds_male, Hess = TRUE)
summary(overall_olr_2_m)

overall_olr_3_m = polr(WLG3_f ~ rand_exp, data = ds_male, Hess = TRUE)
summary(overall_olr_3_m)

summary_table1_m <- coef(summary(overall_olr_1_m))
pval <- pnorm(abs(summary_table1_m[, "t value"]),lower.tail = FALSE)* 2
summary_table1_m <- cbind(summary_table1_m, "p value" = round(pval,3))
summary_table1_m

summary_table2_m <- coef(summary(overall_olr_2_m))
pval <- pnorm(abs(summary_table2_m[, "t value"]),lower.tail = FALSE)* 2
summary_table2_m <- cbind(summary_table2_m, "p value" = round(pval,3))
summary_table2_m

summary_table3_m <- coef(summary(overall_olr_3_m))
pval <- pnorm(abs(summary_table3_m[, "t value"]),lower.tail = FALSE)* 2
summary_table3_m <- cbind(summary_table3_m, "p value" = round(pval,3))
summary_table3_m

#Female

ds_women$WLG1_f <- as.factor(ds_women$WLG1)
ds_women$WLG2_f <- as.factor(ds_women$WLG2)
ds_women$WLG3_f <- as.factor(ds_women$WLG3)

overall_olr_1_f = polr(WLG1_f ~ rand_exp, data = ds_women, Hess = TRUE)
summary(overall_olr_1_f)

overall_olr_2_f = polr(WLG2_f ~ rand_exp, data = ds_women, Hess = TRUE)
summary(overall_olr_2_f)

overall_olr_3_f = polr(WLG3_f ~ rand_exp, data = ds_women, Hess = TRUE)
summary(overall_olr_3_f)

summary_table1_f <- coef(summary(overall_olr_1_f))
pval <- pnorm(abs(summary_table1_f[, "t value"]),lower.tail = FALSE)* 2
summary_table1_f <- cbind(summary_table1_f, "p value" = round(pval,3))
summary_table1_f

summary_table2_f <- coef(summary(overall_olr_2_f))
pval <- pnorm(abs(summary_table2_f[, "t value"]),lower.tail = FALSE)* 2
summary_table2_f <- cbind(summary_table2_f, "p value" = round(pval,3))
summary_table2_f

summary_table3_f <- coef(summary(overall_olr_3_f))
pval <- pnorm(abs(summary_table3_f[, "t value"]),lower.tail = FALSE)* 2
summary_table3_f <- cbind(summary_table3_f, "p value" = round(pval,3))
summary_table3_f

overall_olr_1_m$AIC <- AIC(overall_olr_1_m)
overall_olr_2_m$AIC <- AIC(overall_olr_2_m)
overall_olr_3_m$AIC <- AIC(overall_olr_3_m)
overall_olr_1_f$AIC <- AIC(overall_olr_1_f)
overall_olr_2_f$AIC <- AIC(overall_olr_2_f)
overall_olr_3_f$AIC <- AIC(overall_olr_3_f)

stargazer(overall_olr_1_m, overall_olr_1_f, overall_olr_2_m, overall_olr_2_f, overall_olr_3_m, overall_olr_3_f, 
          type = "html", digits = 2, summary = T, 
          dep.var.labels = c("Lay judge",
                             "Probationer", 
                             "Correspondent"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          keep.stat = c("aic", "rsq", "n"),
          out ="../results/regression_robustness2.htm")

# 7.2 Exclude Failed Manipulation Check ##################################

# 7.2.1 Full Sample =====================================================

ds_treat_s$rand_f <- as.factor(ds_treat_s$rand_exp) 

overall_rev1 <- lm(WLG1 ~ rand_f, data = ds_treat_s)
summary(overall_rev1)

overall_rev2 <- lm(WLG2 ~ rand_f, data = ds_treat_s)
summary(overall_rev2)

overall_rev3 <- lm(WLG3 ~ rand_f, data = ds_treat_s)
summary(overall_rev3)

overall_rev1_d <- glm(WLG1_d ~ rand_f, data = ds_treat_s)
summary(overall_rev1_d)

overall_rev2_d <- glm(WLG2_d ~ rand_f, data = ds_treat_s)
summary(overall_rev2_d)

overall_rev3_d <- glm(WLG3_d ~ rand_f, data = ds_treat_s)
summary(overall_rev3_d)

stargazer(overall_rev1, overall_rev1_d, overall_rev2, overall_rev2_d, overall_rev3, overall_rev3_d, 
          type = "html", digits = 2, summary = T, 
          dep.var.labels = c("Lay judge",
                             "Lay judge (>3)",
                             "Probationer",
                             "Probationer (>3)",
                             "Correspondent",
                             "Correspondent (>3)"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          keep.stat = c("aic", "rsq", "n"),
          out ="../results/regression_manip_s.htm")

# 7.2.2 Sub Samples =================================================

ds_treat_s$gender <- as.numeric(ds_treat_s$gender)

ds_women_s <- ds_treat_s[ which(ds_treat_s$gender == 2), ]

ds_male_s <- ds_treat_s[ which(ds_treat_s$gender == 1), ]

overall_w_rev1 <- lm(WLG1 ~ rand_f, data = ds_women_s)
summary(overall_w_rev1)

overall_w_rev2 <- lm(WLG2 ~ rand_f, data = ds_women_s)
summary(overall_w_rev2)

overall_w_rev3 <- lm(WLG3 ~ rand_f, data = ds_women_s)
summary(overall_w_rev3)

overall_m_rev1 <- lm(WLG1 ~ rand_f, data = ds_male_s)
summary(overall_m_rev1)

overall_m_rev2 <- lm(WLG2 ~ rand_f, data = ds_male_s)
summary(overall_m_rev2)

overall_m_rev3 <- lm(WLG3 ~ rand_f, data = ds_male_s)
summary(overall_m_rev3)

stargazer(overall_w_rev1, overall_m_rev1, overall_w_rev2, overall_m_rev2, overall_w_rev3, overall_m_rev3, 
          type = "html", digits = 2, summary = F, 
          dep.var.labels = c("Lay judge",
                             "Probationer",
                             "Correspondent"),
          covariate.labels = c("Mixed treatment",
                               "All female treatment"),
          omit.stat = "F", star.cutoffs = c(0.05, 0.01, 0.001),
          out ="../results/regression_manip_s_gender.htm")

# 8. End #################################################################